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Competitive birth-death processes often exhibit an oscillatory behavior. We investigate a par- 
ticular case where the oscillation cycles are marginally stable on the mean-field level. An iconic 
example of such a system is the Lotka-Volterra model of predator-prey competition. Fluctuation 
effects due to discreteness of the populations destroy the mean-field stability and eventually drive 
the system toward extinction of one or both species. We show that the corresponding extinction 
time scales as a certain power-law of the population sizes. This behavior should be contrasted with 
the extinction of models stable in the mean-field approximation. In the latter case the extinction 
time scales exponentially with size. 
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I. INTRODUCTION 

Understanding stochastic population dynamics is an 
important pursuit in the biological sciences [l], [H, 0, 0] • 
Mathematical modeling of such dynamics allows for bet- 
ter understanding of biodiversity and species extinction. 
Such modeling becomes especially important because of- 
ten the relevant time scales make direct measurements 
difficult. One of the most basic relationships that can be 
used to study such dynamics is the predator-prey rela- 
tion. In such a system, one species reproduces by killing 
the other. An individual of the prey species replicates 
at a constant rate. Individual predators die at a con- 
stant rate and replicate only at the expense of the prey. 
Although the most obvious application of such a system 
is two organisms, the predator-prey relation can also be 
used to study other systems. 

The original work by Lotka [|[ and Volterra @ showed 
that such a system results in oscillations of both pop- 
ulations. Stochastic simulations can be used to better 
understand such a system. In a system without spatial 
degrees of freedom, the Lotka-Volterra interaction invari- 
ably results in an extinction event in which either the 
predator species or prey species goes extinct. This de- 
parture from the original results can be understood as a 
result of stochasticity associated with the discreteness of 
the populations. Such behavior has been observed in, for 
example, the cyclic Lotka-Volterra system Q. Under- 
standing this departure from mean-field dynamics pro- 
vides a challenging problem in non-equilibrium statistical 
mechanics. 

The unique feature of the Lotka-Volterra model is the 
presence of an "accidental" first integral of the mean-field 
equations of motion. As a result, all mean-field trajecto- 
ries evolve on closed orbits. These type of dynamics are 
marginally stable, since fluctuations in any direction are 
neither damped nor amplified. Such fluctuations origi- 
nate from intrinsic demographic stochasticity along with 
the discreteness of the populations; they lead to a slow 
diffusion between the mean-field orbits. Even large devi- 
ations from mean-field expectations, such as extinction, 
may be viewed as the accumulation of many small step 



fluctuations in the radial direction. This should be con- 
trasted with reaction systems that have a stable fixed 
point or limiting cycle. In those system the large devia- 
tions proceed only along very special instanton paths in 
the phase space [1, Q. Due to this difference, the ex- 
tinction time in marginally stable systems exhibit power 
law dependence on the two populations sizes, instead of 
being exponentially long as in the case of (meta)stable 
models. 

This work gives a formulation of the problem using 
the Fokker-Planck equation. Using time scale separa- 
tion between fast angular and slow radial motion, the 
inherently two-dimensional problem is reduced to a one- 
dimensional one. The latter is the problem of diffusion 
with a specific radius dependent drift. We then solve the 
first passage problem for this effective ID problem and 
characterize the extinction probability in the long and 
short time limits. We rely on extensive comparison of 
the analytic results with the stochastic simulations. We 
achieve a quantitative agreement between the two, which 
is in all cases is better than 5% and may reach an accu- 
racy of 0.5%. 

Our main result may be formulated as follows: for 
generic parameters and initial conditions the typical 
number of cycles C the system undergoes before going 
extinct scales as 

C cx A s 3 / 2 x N~ 1/2 , 

where > N s are the sizes of the dominant and sub- 
dominant populations correspondingly. This result im- 
plies a number of surprising consequences, which were all 
confirmed in simulations. For example, it predicts that a 
further increase of an already dominant population only 
accelerates the total extinction. It also shows that some 
very different systems behave virtually indistinguishably 
vis-a-vis extinction, if their C-numbers are the same. For 
the symmetric scenario Nd — N s — N, we find C cx N in 
agreement with Ref. 

The outline of this paper is as follows: in section ITT1 we 
present the mean-field dynamics of the Lotka-Volterra 
system. Section Mil presents some of the results of ex- 
tensive Monte Carlo simulations. An analytic approach 
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to understanding the problem is presented in section [TV] 
Finally, the results are discussed in section fVl 

II. MEAN-FIELD THEORY 

In a basic predator-prey system, there are two popula- 
tions. The predator species has a death rate, a, and the 
prey has a birth rate, In addition, there is a cross re- 
action where a predator consumes a member of the prey 
population in order to reproduce. This occurs at rate A. 
The reaction scheme can be summarized as follows: 

F'O; R^2R; F + R ^ 2F , (1) 

where F signifies a predator ("fox") and R signifies a 
prey individual ("rabbit"). 

In the mean-field approximation one neglects the dis- 
creteness of the populations and models the system with 
deterministic rate equations. If q± and 52 are taken to 
be continuous variables representing the predator and 
prey populations, the dynamics of these two variables 
are given by the following equations 

4i = -om + , 

q-2 = Mi - Agi<?2 • (2) 

The rate of change of the predator population contains a 
death term proportional to the predator population and 
a birth term proportional to the size of both populations. 
Likewise, the rate of change of the prey population has 
a birth rate proportional to the prey population and a 
death term proportional to both. Some features of these 
dynamics arc immediately evident. There are three fixed 
points. These correspond to (ffrjffe) = (0,0), (0, 00), and 
(/x/A,(j/A). The first point corresponds to the trivial case 
of extinction of both species. The second fixed point is 
the result of predator extinction and the prey population 
growing exponentially. The third is the coexistence fixed 
point, where the stable populations of the predator and 
prey are Ni = ///A and N2 = er/A. 

For a given initial condition, the populations evolve 
along a closed orbit in predator-prey space. The orbits 
are closed due to the existence of an " accidental" integral 
of motion in the mean-field equations of motion ([2]) 

G = Agi — fi — fihi (qiX/fi) + Aq 2 — a — a In (q 2 X/a). (3) 

The definition of G is chosen such that G = corre- 
sponds to the coexistence fixed point, while G — > 00 cor- 
responds to large amplitude cycles closely approaching 
the two axes. Figure [1] shows orbits for various values of 
the integral of motion, G, for the case Ni = N2 — 100. 
The presence of the integral of motion makes all cycles 
marginally stable. Indeed, a small fluctuation may shift 
the system from one orbit to a neighboring one. Since the 
new orbit is also a stable solution of the mean-field equa- 
tions of motion, there is neither a restoring force, trying 
to compensate for the fluctuation, nor amplification of 
the fluctuation. 
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FIG. 1: Orbits of constant G = (0.01, 0.1, 0.4, 1, 1.7, 2.7, 4.2) 
in units of ^/cr/i. The evolution proceeds clockwise around the 
mean-field fixed point of Ni — N2 = 100. 
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FIG. 2: Typical run of the stochastic simulation of the model 
HI for N = 100 and e = 1. 



At small G, the mean-field orbits are ellipses. The fre- 
quency of these elliptical orbits approaches a constant 
value, l/y/ajl, as G approaches zero. This provides a 
natural time scale for the problem. By rescaling time to 
be measured in these units, it is possible to reduce the 
number of parameters in the problem from the original 
three reaction rates to two parameters, which are conve- 
nient to choose as 
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Here N represents an effective system size, while e rep- 
resents the asymmetry between the predator and prey 
populations. Throughout this paper we shall be inter- 
ested in the limit of large system size N ^> 1. By the 
reasons explained below the asymmetry parameter e will 
be restricted to the interval iV~ 1/2 < e < TV 1 / 2 . 



III. STOCHASTIC SIMULATIONS 

The mean-field approximation fails to accurately por- 
tray the actual evolution of the reaction system, Eq. (JTJ) . 
As previously mentioned, the mean-field solution does 
not take into account the stochasticity associated with 
individual birth-death events or the discreteness of the 
populations. Results from Monte Carlo simulations of 
this reaction scheme demonstrate the failure of the mean- 
field. 

Stochastic simulations were done as follows. The ini- 
tial populations were taken to be at the coexistence fixed 
point. Time was discretized into small steps of size St. 
The time step, St, was chosen so that the probability of 
having any change in population size during St was small, 
i.e. St <C \jNJW\i. For each time step, the number of 
prey births and predator deaths was calculated randomly 
from a binomial distribution with success rates fiSt and 
oSt respectively. The number of consumption reaction 
events was calculated from a binomial distribution based 
on XSt and the number of predator/prey pairs. This was 
repeated until one of the populations went extinct. Fig- 
ure shows an example of such a simulation. As in the 
mean-field case, the system rotates clockwise about the 
coexistence fixed point. As time progresses, however, the 
system unwinds from the fixed point, eventually hitting 
either the q\ or q2 axis. From there, the system rapidly 
progresses toward one of the extinction fixed points. For 
a typical simulation, the system rotates around the fixed 
point many times before going extinct. 

Since such a simulation invariably ends in the extinc- 
tion of one or both species, it is interesting to analyze the 
chance of the system being dead as a function of time. 
For a given set of initial conditions, it is possible to de- 
termine this extinction probability by repeatedly running 
a stochastic trial. Figure [3] shows the result of 100,000 
stochastic simulations using the conditions of the simula- 
tion presented in Fig. [5J As could be expected, at short 
time scales there is very little chance for the system to 
be extinct. As t grows the probability of extinction ap- 
proaches unity. At long time scales, the convenient quan- 
tity for calculation is not the extinction probability, but 
the survival probability, this being the likelihood of the 
system still being alive at a given time, t. The logarithm 
of the survival probability appears to be linear in time, 
Fig. 2J As a result, the survival probability is decaying 
exponentially with a characteristic time 77, 

P surv (t) = 1 - P ext (t) oc e-^ Tl . (5) 

Figure O shows the dependence of 77 on N. One observes 



H 0.8 

-rH 



To 0.6- 
u 

Cm . 

g 0.4 - 

-H 
-p 
O 

■h . 2 - 

+J * 
X • 
w • 

50 100 150 200 

t 

FIG. 3: Extinction probability in time t from 10 5 simulation 
trials (N = 100, e = 1). Time is in units of l/y/ajl. 
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FIG. 4: Logarithm of the survival probability at long times 
for the data presented in Fig. [3] 

a linear growth of the characteristic time 77 with increas- 
ing N at N 3> 1. This agrees with the results observed 
by Reichenbach, et al. Q for the cyclic Lotka-Volterra 
system. This linear dependence suggests the following 
representation for Ti(N, e) in the limit N ^> 1: 

where the rescaled extinction rate, Eq, depends only on 
the asymmetry, e, but not on N. The fit of Fig. [5] gives 
an observed value of E (l) = 2.05. 

We focus now on the role of the asymmetry parameter 
e. Figure [5] plots the extinction probabilities versus time 
for e = 2 and e = 1/2. The plots show virtually identical 
behavior. In particular, the similarity in the long time 
decay suggests a symmetry in Eq between e and 1/e. Fig- 
ure [7] shows a plot of the observed Eq vs. the logarithm 
of e in stochastic simulation, confirming that 

E (e) = E (l/e). (7) 

The minimum of Eq corresponds to e = 1. Away from 
this point one observes E (e) — 0.97(max{e, 1/e}) . 

Unlike the long times, the short time behavior is not 
universal and depends on the choice of the initial condi- 
tions. For the initial populations chosen close to the coex- 
istence fixed point (N±,N2) it is exceedingly unlikely for 
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FIG. 5: Time constant t; of exponential decay, Eq. ([5]), versus 
N for e = 1. 
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FIG. 8: Logarithm of the extinction probability at short times 
for the data of Fig. [3] 
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FIG. 6: Extinction probabilities for e = 2 and e = 1/2; 
N=100. 



linearly with AT for N ^> 1. It is thus convenient to pa- 
rameterize t s (N, e) by an iV-independent quantity Xo(e) 
as 



r s (AT,e) 



iV 



(9) 



This parameterizations of r s puts Eq.® in a form that 
is reminiscent of a standard diffusion propagator. As was 
the case for Eq(e), Xo(e) also shows symmetry between e 
and 1/e, i.e. X (e) = X (l/e). From the simulations one 
observes Xq(1) = 2.09 at the maximum. 



IV. ANALYTIC APPROACH 



the system to drift all the way to extinction in a short 
time interval. Plotting the logarithm of the extinction 
probability vs. inverse time shows a linear dependence 
at small times. This can be seen for e = 1 in Fig. [5J Ex- 
tinction probability in small times is thus exponentially 
small and appears to have a functional form of 

P ext cx e" T */*. (8) 

As in the long time limit, the time constant t s (N, e) scales 
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FIG. 7: Plot of E vs. In e ; N = 100. 



A. Master and Fokker-Planck Equations 

The full behavior of the reaction model {J) can be ana- 
lyzed by employing a probability distribution and study- 
ing its dynamics. Define such a probability distribu- 
tion P(m, n; t) as the probability of the system having m 
predators and n prey at time t, where m and n are both 
non-negative integers. This yields the following master 
equation for the reaction scheme of Eq. |T]) 

dtP(m, n; t) = a[(m + \)P(m + 1, n) — mP(m, n)\ 
+ H[(n-l)P(m,n-l)-nP(m,n)] (10) 
+ A[(m — l)(n + l)P(m — 1, n + 1) — mnP(m, n)]. 

This equation can be rewritten using integer shift oper- 
ators, defined as 



E, 



k.l 



kd m +W n 



(11) 



to give 



dtP{m,n]t) = a(Exfi — l)mP(m, n; t) 

+ n{E Q - 1 -l)nP{m,n-t) (12) 
+ A(-ELi,i — l)mnP(m, n; t) . 

An important distinction of the models with 
marginally stable cycles is that a large deviation (such 
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as an extinction) may proceed in a sequence of small 
steps. A small fluctuation leads to a mean-field like 
evolution along a new stable orbit until another small 
fluctuation shifts the orbits again, etc. As a result, a 
path to extinction is a random diffusion in population 
space. This should be contrasted with models with sta- 
ble limiting cycles or an attracting fixed point [l(| EH , 
where small fluctuations do not accumulate and extinc- 
tion proceeds only along a very specific (instanton) tra- 
jectory. On the technical level this observation implies 
that the gradients d m>n may be considered as small 
~ 1/N (this is usually not the case on the instanton tra- 
jectory H, !, E3, 0, 0, GH) and thus the shift operators 
(fTTjl may be expanded up to the second order. This pro- 
cedure leads to the Fokker-Planck (FP) equation, which 
in the present context is justified by the Van-Kampen 
expansion over the system size N [161] . Proceeding this 
way, one finds 



d t P = o 
+A 



5* 



8a 



1 



qiP + fi 
1 



-d, 



-d q2 - d qi + -d qi -d qi d, 



</2 



q 2 P 
qiq 2 P- (13) 



This equation along with Eq. ([2]) suggest a change of vari- 
able such that the new variables Qi ~ Inqi. This is ac- 
complished through the following transformation: 



A 



92 



A 



(14) 



These variables present some advantage over the initial 
ones. Extinction events now occur at Q\ = — oo or 
Q2 = —00 instead of at q\ = or <?2 = 0. The coex- 
istence fixed point has been moved to the origin. As part 
of this transformation, time is rescaled into the prob- 
lem's natural units, l/y/ajl. The FP equation no longer 
depends on the three reaction rates; it depends only on 
N and e. In the new coordinates the mean-field integral 
of motion takes the form 

G = -(e^ - 1) - Qi + e(e Q ^ £ - 1) - Q 2 . (15) 
e 

It provides a natural radial coordinate. The coexistence 
fixed point is G = 0, while extinction corresponds to 
G = 00. Figure [5] shows mean-field orbits in the trans- 
formed coordinate system. Larger orbits correspond to 
larger values of G. The most essential advantage of the 
new variables is that the mean-field equations (J2]) ac- 
quire the Hamiltonian structure, where Q± and Q2 form 
a canonical pair 



h = -l + e^l' =d Q2 G; 
h = l-e^=-d Ql G. 



(16) 



Since G(Qi,Q2) serves as the Hamiltonian, it is mani- 
festly conserved on the solutions of the mean-field equa- 
tions of motions. 



The probability distribution is transformed in the new 
coordinate system so as to include the Jacobian of the 
transformation 



W(Qi,Q 2 ;t) = q 1 q 2 P(qi,q2;t), 



(17) 



where qi and q 2 are substituted from Eq. (fl"4"|) . In the new 
coordinates, the Fokker-Plank equation (fT3")) becomes 

d t W — —V • J, (18) 
where the divergence is defined as 

v = (d Ql ,d Q2 ). (19) 

The probability current in Eq. (|18|) consists of two parts 
J = J MF + J D . (20) 

The mean-field motion along the orbits of constant G 
is due to J MF , see Eq. (fTo]) , while the radial diffusion 
between the orbits is due to J D . The mean- field current 
is given by 

J* F = (-1 + e Q2 / e )W = (d Q2 G)W ; (21) 



J. 



MF 



)W 



\d Ql G)W. 



(22) 



The diffusive current is found from Eq. (|13p as 



D 



1 

~2N 



2N 



(e 



-eQl 



<?2 



eQ 



')d Ql W-d Q2 W ■ (23) 



(e — + e c 



_ Q-i 



)d Q2 W-d Ql W . (24) 



The diffusive current is suppressed by a factor of N rel- 
ative to the mean-field current. Provided the system is 
sufficiently large, i.e. AT 3> 1, this means that the angular 
motion should be much faster than the radial motion. 



B. Reduction to One Dimension 

As mentioned, the mean-field constant, G, provides 
a natural radial coordinate. The corresponding angular 
coordinate evolves far faster than the radial one. This 
time-scale separation represents an opportunity to turn 
this two-dimensional problem into a one-dimensional one. 
The method used has been successfully employed in the 
analysis of spin- torque switching (l7| . Since Q\ and Q2 
form a canonical pair on the mean-field level, it is possi- 
ble to transform them into action-angle variables (/, a) 
where the action is an integral of the mean-field motion, 
i.e. G = G(I). It is more convenient to use G itself, 
rather than the canonical action, /, as the radial variable. 
One should be aware though that the change of variables 
(Qi,Q2) — * (G,a) involves a Jacobian, discussed below. 
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FIG. 9: Orbits of constant G in the new coordinates 



We shall assume now that due to the fast angular mo- 
tion the probability distribution W(G, a; t) rapidly equi- 
librates in the angular direction and at the long time 
scales depends on the radial variable only, i.e. W — 
W(G;t). Under this assumption it is possible to elimi- 
nate the angular dependence from Eq. (fl~8|) by integrating 
the continuity equation (fTS)) over the area of a mean-field 
orbit with a fixed G 



d t W(G)dQ 1 dQ 2 



V ■ JdQ 1 dQ 2 (25) 



We apply the divergence theorem to the right hand side 
and change the coordinates of integration on the left hand 
side. The unit vector, h, is normal to the line of constant 
G, and dl is an infintismal length along the G orbit 



dtW(G) 



G 



dQx dQ 2 8Q 2 dQi 



dG da dG da 



G 



dGda =—<t> J ■ ndl 
(26) 

The Jacobian can be reduced to \dt/da\ due to the Hamil- 
tonian equations on Qi and Q 2 , see Eq. (|16p . leading to 



d t W{G) 



dt 
da 



dadG 



J ■ hdl. 



(27) 



G 



Integration over a gives the period of the mean-field rev- 
olution around the orbit T(G) 



T(G) d t W{G)dG 



J ■ ndl . 



(28) 



Finally, differentiation with respect to G yields the radial 
FP equation 



T(G)d t W{G) = ^ 



- - </> J ■ ndl 



(29) 



We now wish to understand the integral of the current 
in this equation. The mean-field portion of the current, 
J , is perpendicular to n. It therefore makes no con- 
tribution to this integral 



J 



'MF 



0. 



(30) 



This leaves integration over the diffusive current J D , 
which is first-order in derivatives and proportional to 
1/N. Since we are assuming that W is independent of a, 
this implies that the diffusive current should be propor- 
tional to daW . The integral along the orbit may then be 
written in the form 



J°.ndl = ±DlGr— 



dW 
~d~G 



(31) 



which is in essence the definition of the effective diffusion 
parameter, D(G). The resulting FP equation takes the 
form 



T(G)d t W(G) = ^ 



1 D(G) dW{G) 
N D{G) dG 



(32) 



where the two functions D(G) and T(G) may be evalu- 
ated for any mean-field orbit G. They both are indepen- 
dent of N, but do depend on e. We evaluate both these 
quantities in the Appendix. At small G <C 1, the period 
is T = 2-7T (which was our initial motivation of choosing 
these units of time) and D(G) oc G. Introducing variable 
R = VG, so 8g = (l/2R)du, one reduces the right hand 
side of Eq. ([32]) to be ~ R- 1 d R (Rd R W), which is the ra- 
dial part of the standard 2d diffusion equation. At large 
G the period grows linearly, T(G) ~ G, while the diffu- 
sivity D(G) grows exponentially. The latter fact is due 
to the two sharp maxima in the current J D which take 
place around the two "arms" of the mean- field orbits, 
Fig. [5] along the negative directions of the two axis. Both 
T(G) and D{G) are symmetric with respect to e — > 1/e, 
rendering the corresponding symmetry to all the results 
obtained upon averaging over the angular motion. 

The boundary conditions for Eq. (|32p are as follows: 
since the system is incapable of moving from the extinct 
state to a live one, there is an absorbing boundary con- 
dition as G — * oo. Since G is a radial coordinate, the 
current must disappear at G = 0. This gives the follow- 
ing boundary conditions 



lim W(G) = 0: 

G^OG 



G 



dW(G) 



dG 



G=0 



= 0, 



(33) 



where we have employed D(G) ~ G at G — > 0. At 
large G, the diffusivity D(G) grows exponentially, see 
Appendix. This allows the system to diffuse to G = oo in 
finite time which suggests that the spectrum of Eq. (|52"]) 
with the boundary conditions (j33|) is discrete, despite 
the equation being formulated on the infinite interval 
G G [0, oo). To see this fact most clearly it is useful 
to introduce one more change of variables. 
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C. Reduction to the finite interval 

To motivate the new variable, let us perform a semi- 
classical analysis of Eq. (f32|) . To this end we represent 
the probability distribution as 



W(G;t) oce- NS ^ 



(34) 



and assume for the moment that S ~ 0(1) (this assump- 
tion is indeed true at short time scales, but breaks down 
at t ~ N). Then to the leading order in N, Eq. ([32]) 
becomes 



T(G)d t S(G;t) = D(G) 



dS{G;t) 
dG 



(35) 



which may be viewed as the Hamilton-Jacobi equation 
with the Hamiltonian 



H(G,P G ) 



D(G) 
T{G) 



(36) 



It is convenient to make a canonical transformation from 
(G,Pq) — > (X,Px) such that the Hamiltonian in the 
new coordinates takes the form H(X, Px) = P x - This is 
accomplished by 
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FIG. 10: X(G) for e = 1. At G — > oo the function converges 
to Xo = 2.39 
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(37) 



Unlike G, the new radial variable X is bounded. Indeed, 
due to the exponential growth of D(G) at large G the 
integral in Eq. (|3"T|) converges. For the case e = 1, X(G) 
is plotted in Fig. [Till exhibiting convergence to X — 2.39. 

Although we motivated the change of variables G — > 
X(G) by the semiclassical analysis of Eq. (|32p . one may 
go back to the full FP equation (|3"!?|) and perform the 
variable change ([37]) exactly. The result is 



d t W= — 



V J D{X)T{X) ^ {^ D ^ T W 



(38) 

This equation is defined on the finite interval X 6 [0, Xo]. 
The boundary conditions ([33]) take the form 



W(X ;t) =0; 



X 



dW(X;t) 



dX 



x=o 



0, 



(39) 



where we take into account that y/ D{X)T(X) ~ X at 
X — * 0. Equation (|38|) has framed the problem so that 
it no longer depends on T and D separately, but rather 
o nly on \JTD as well as the constant Xq. A plot of 
y/T(X)D(X) is shown in Fig. [JTJ 



FIG. 11: Numerically calculated V / T(X)D(X) for e = 1. The 
function diverges at X — Xo = 2.39 



D. Long time dependence of the extinction 
probability 

The long time behavior of the system can be analyzed 
via an eigenvalue problem on Eq. (|38p with the bound- 
ary conditions (|39p . Since Xq is finite, its spectrum is 
discrete. We shall look thus for a solution of the FP 
equation in the form 



a n w n (X) i 



-E n t/N 



(40) 



where a n are constants depending on initial conditions 
and w n (X) and E n are solutions of the following eigen- 
value problem 



Hw n (X) =E n w n {X). 

Here the operator H is defined as 
-1 d 



H 



X / D(X)T(X) dX\^ D{X)T{X) dX 



(41) 



(42) 



The N dependence has been explicitly removed from the 
eigenvalues. The only remaining dependence of E n is 
on the asymmetry e. The survival probability is given 
by P S urv{t) — J Q X ° dX W (X ; t) . At long time scales, the 
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FIG. 12: Eo numerically calculated via matrix diagonalization FIG. 13: Function X 1 (e). 

vs. number of rows in matrix 



only contributing eigenstate is the one with the smallest 
eigenvalue, Eq, and thus the survival probability decays 
with the characteristic time scale r; = N/Eq(c), 



J*surv (0 CC e 



(43) 



This is what was observed in Fig. [5J From the fit of this 
figure, the observed value from stochastic simulation is 
E Q (1) = 2.05. 

In order to compare it with our analytical approach 
we discretize the interval [0, Xq] and represent the linear 
operator H, Eq. (|42| . as a matrix. We then diagonalizc 
it numerically, using expressions for T and D from the 
Appendix. The lowest eigenvalue, Eq, converges quite 
rapidly with decreasing discretization step of X; 100 rows 
is sufficient to calculate Eo to within 1% of the convergent 
value. Figure [T^] shows the convergence of E as the 
number of matrix rows is increased. This procedure gives 
for the operator H a lowest eigenvalue of Eq{1) — 1.95, 
which agree with the Monte Carlo simulations within 5% 
accuracy. The discrepancy can be reduced even father 
by taking into account the finite size effect. For finite 
values of N, it is not necessary to diffuse all the way to 
G = co, but only to a value of G which corresponds to a 
single remaining individual at a minimum of one of the 
populations. At this point the fluctuations will drive the 
system to extinction with probability close to one. Such 
a cutoff G ext may be estimated, using Eq. ©, as 



e- 1 (ln(Ar/e)-l); e>l 
e(ln(iVe) - 1) ; e < 1 



(44) 



For our simulations with N — 100 and e = 1 this gives 
Gext = 3.62. Integrating Eq. ([3T| only up to G ext gives 
Xq = 2.08 (instead of Xq — 2.39 for the infinite interval). 
Using this truncated Xq as the cut-off for the matrix di- 
agonalization of H gives Eq(1) = 2.06, within 0.5% of 
stochastic simulations. Since G ex t depends on the sys- 
tem size N only logarithmically, it is computationally 
unfeasible to eliminate the finite size effect in stochastic 
modeling. 

We turn now to the e dependence away from e = 1. 
Let us discuss the e > 1 case, i.e. N 2 > N% prey domi- 
nated system (the predator dominated scenario may be 



analyzed in the same way). In this case it is almost cer- 
tain that predators go extinct first. This is because the 
diffusive current toward predator extinction, Eq. (|A.7[) . 
is exponentially bigger than that toward prey extinction. 
Neglecting the latter, one observes from Eq. (|A.7p that 



D(eG) 



This implies that the integration 



D 

interval contributing to Xq, Eq. (|3"T)) . is effectively lim- 
ited to < G < 1/e < 1. In this interval the pe- 
riod T(G) « const. Rescaling variables in Eq. (|37f as 
eG -> G, one finds that X (e) - 1/e for e > 1, Fig. [ED 
Correspondingly, D(X) — D(eX) and after rescaling 
(X -> X in Eq. (@2]) one observes that E n (e) ~ e 2 . Fi- 
nally, one finds for the characteristic extinction time of 
an asymmetric model 



N 



Eo(e) 



1.03 iV(max{e, 1/e})" 



(45) 



where the numerical factor is obtained through numer- 
ical diagonalization of the H operator. Figure [TJ] plots 
the observed values from Monte-Carlo simulation fit with 
our analytic prediction. Since the approach relies on the 
separation of time scales between the fast angular mo- 
tion and the slow radial one, it requires t\ > 1. This 
leads to the restriction on the asymmetry parameter: 
AT -1 / 2 < e < TV 1 / 2 , stated in section OH Outside of this 
interval it takes about a period of one small revolution 
for the system to go extinct. 



E. Short time dependence of the extinction 
probability 



In the short time limit, the semi-classical analysis 
presented in subsection IIV CI should be accurate. In- 
deed, extinction in time t <C t s is an exponentially 
rare event thats probability is convenient to represent as 
W = e~ NS . Here the action S(G;t) is a solution of the 



Hamilton- Jacobi equation (|35p with the initial condition 
5(0; 0) = and G — > oo at time t . After the canoni- 
cal transformation (G,Pa) — > (X,Px), the Hamiltonian 
acquires the form H(X, Px) — P\ and the classical equa- 
tion of motion is X = 2Px ■ We need its solution reaching 
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FIG. 14: Function Eo(e); dots - results of stochastic simula- 
tion, full line - operator H diagonalization. 



FIG. 15: Extinction probability of the two models: crosses 
N = 100, e = 5; triangles N = 400, e = 10. In both cases 
ti » 4.0, while r s « 8.8. 



X = at time i. The corresponding action is 



S(X ;t) 



(46) 



resulting in the following form for the short time scale 
behavior of the extinction probability 



P. 



ext ^ ^ 



-JVX 2 /(4t) 



(47) 



This is exactly what was observed in stochastic simula- 
tions. The fit from Fig.[4]gives Xo = 2.09 for e = 1, while 
evaluating X from Eq. (|37|) results in A = 2.39. Again, 
the majority of the difference between these two values 
can be eliminated by using the value Xq that is corrected 
for the finiteness of TV, cf. Eq. (HH). This was calculated 
in the previous section to be Xq = 2.08, in much better 
agreement with the simulations. The e dependence of the 
short time scale r s = AAq(£)/4 follows from the depen- 
dence X (e) discussed in the previous section. Thus, one 
finds that away from the symmetric point e = 1 



t s = 2.2 A(max{e, l/e})~ 



(48) 



where the numerical constant is obtained through numer- 
ical integration of Eq. ([37)) . 

One may argue that at the very smallest time scales 
the reduction of the initial Master equation (fTU|) to the 
FP equation JTjJ) may not be justified. As a result at such 
small times either Eq. (|47|) . or Eq. (|48|) may be violated. 
Although this is a potentially valid concern, we were not 
able to go to sufficiently short times (or sufficiently large 
N) to detect any sizable deviations of Monte Carlo results 
from analytical predictions, Eqs. (|47|) , ([48]) . 



V. DISCUSSION 

We have investigated extinction due to intrinsic 
stochasticity in the Lotka-Volterra model (JTJ) . To this end 
we have introduced two characteristic times: (i) the uni- 
versal scale 77, which characterizes exponential decay of 
survival probability at long times; (ii) the non-universal 



scale t s specific to the choice of initial condition close 
to the coexistence fixed point, which characterizes rise 
of the extinction probability at short times. Since both 
these scales depend on the system parameters in exactly 
the same way and differ from each other only by a fac- 
tor close to two, t s w 2.2t;, we shall restrict ourselves to 
discussions of the time 77, which is independent on the 
choice of the initial conditions. All our results are valid 
in the asymptotic limit of large system size Ni t 2 S> 1. 

We consider first the asymmetric case. Recalling the 
definition of the parameters, Eq. Q, and employing 
Eq. ([45]) one finds 



N. 



3/2 



n = 



N. 



1/2 



(49) 



where the size of the dominant population is = 
max{iVi, N2}, the size of the subdominant one is N s = 
minjiVi, N2}, and time is measured in the natural units, 
which is the inverse frequency of the small cycles 1/ \faji. 
This is a remarkable scaling relation, which predicts e.g. 
that the extinction time is shortened with increasing the 
size of the dominant population. Countcrintuitively, in- 
creasing abundance of dominant "rabbits" accelerates the 
extinction of subdominant "foxes"! To check this pre- 
diction we performed stochastic modeling of two prey- 
dominated models, which according to the scaling of 
Eq. (|49| ought to go extinct in the same relative time. 
The results are presented in Fig. [T5] Since our method 
involved assumption that the angular motion is faster 
than the radial one, Eq. (|49|) may be trusted as long as 

Outside of this interval 



1/3 



77 > 1, i.e. N d > N s > AT 
of the parameters the extinction time is about one (in 
relative units). 

Returning to the absolute scale of time and recalling 
that N± = /i/A, while A^ = cr/A, one may rewrite the 
extinction time, Eq. (|4"9"1) , as 



= m/OA) ; <r > V: 
" s cr/ipX) ; (i><7. 



(50) 
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The first line here is the prey-dominated case, while the 
second line the predator-dominated one. Again some- 
what counterintuitively, increasing "rabbits" birth-rate 
accelerates their extinction in the "fox" -dominated word. 
In the symmetric case N d = N s = N ^> 1, we find 



10 



n = 0.51 N 



0.51 

A 



(51) 



where the first result is in the relative time scale, while 
the second in the absolute one. The linear scaling of the 
nearly symmetric model with the system size is in agree- 
ment with the results of Ref. [7], obtained for a closely 
related cyclic model. The factor close to a half in compar- 
ison with the asymmetric case, Eq. (|49|). admits a simple 
interpretation. In the asymmetric case the diffusive cur- 
rent toward the extinction of the dominant population is 
exponentially smaller than that toward the extinction of 
subdominant species and may be neglected. In the sym- 
metric case the two currents are exactly the same, making 
the extinction time twice shorter. How close to the sym- 
metric point the system has to be for Eq. (|5Tj) to hold? 
Using Eq. (|A.7|) and taking characteristic value of G from 
Eq. (|44[) one may estimate the corresponding interval of 
parameters as |e-l| < 1/lniV, i.e. N d — N s < N d /\nN d . 
This means that in the limit of large populations the sym- 
metry condition is rather restrictive and a generic system 
most likely obeys the asymmetric scaling. 

The natural extension of our study is inclusion of spa- 
tial degrees of freedom. The spatial extension of the sys- 
tem is capable of stabilizing the system and increasing 
the extinction time 18]. Even in a 2-site system, extinc- 
tion time can be substantially longer than in the zero- 
dimensional case presented here [19(. Understanding of 
such a stabilization mechanism is crucial for an accurate 
description of the phase transition between the absorbing 
extinct phase and active coexistence phase, exhibited by 
the model on a thermodynamically large d-dimensional 
lattice [H. 

We are indebted to M. Dykman and B. Meerson for 
numerous illuminating discussions. This research is sup- 
ported by NSF Grants DMR-0405212 and DMR-0804266. 



APPENDIX: EVALUATION OF D(G) AND T(G) 

In the limit G <C min{e, 1/e}, an orbit of constant G 
is an ellipse. Both parameters, D(G) and T(G), may be 
found exactly in this case 

D{G) = 2vrG(e + l/e); T(G) = 2tt . (A.l) 

Equation (152")) takes the form 



d t W 



+ l/e a 
N dG 



G 



dW 
~dG 



(A.2) 



Changing variables as G = R 2 near the mean-field fixed 
point gives the radial part of the two-dimensional diffu- 




FIG. 16: Numerically calculated T(G) for e = 1 fit with ana- 
lytically predicted T{G) 
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FIG. 17: Numerically calculated D(G) for e = 1 fit with 
analytically predicted D(G) 



sion equation with diffusion constant of (e + 1 /e) /AN 

nTT , e+l/e 1 d ( dW\ , k . 

d * w = ^iiaR{ R lm)- (A - 3) 

The large G limit can also be estimated. The diffu- 
sive current J D ■ h has two maxima corresponding to the 
minima in one of the two species populations. These 
maxima are located at Q2 = 0, Q\ sa — G — 1/e and 
Qi = 0, Q2 ~ — G ~ e. Expanding near these two points 
the currents (f2"5|) and (|2"4"|) and evaluating the integral in 
Eq. (J3TJ) one finds 



/).f;, = j-l e + (i + I) 



\ e eG 



+ J\ (e+(l + e 2 ) 1/2+1/e2 )e GA - (A.4) 

The majority of the orbital period is spent in the third 
quadrant. In this quadrant, Q\ ss —1 and Q\ varies from 
— G to 0. This gives for the orbital period 



T(G) = G . 



(A.5) 



For the purposes of numerical diagonalization of the 
H operator we use the following interpolating function 
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accurate in both the large and small G limits /tt ^ e _|_ (i _|_ £ 2-ji/2+i/c 2 ^ ( e G/e _ gy e _ 

T(G) = 2tt + G. (A.6) V " 

D{G) — 2ttG(€ + 1/e) (A. 7) Figure H7I shows the numerically calculated values for 

/— / -, \ D(G) fit with this interpolated D(G). Figure [TBI shows 

Me+(l + ^) 1/2+£2 J (e eG -eG-l) the same for T(G). 



[1] M. Bartlett, Stochastic Population Models in Ecology and 
Epidemiology (Wiley, New York, 1961). 

[2] H. Andersson and T. Britton, Stochastic Epidemic Mod- 
els and Their Statistical Analysis (Springer, New York, 
2002). 

[3] O. Diekmann and J. Heesterbeek, Mathematical Epi- 
demiology of Infectious Diseases: Model Building, Anal- 
ysis, and Interpretation (Wiley, Chichester, 2000). 

[4] D. Daley and J. Gani, Epidemic Modelling: An Introduc- 
tion (Cambridge University Press, Cambridge, 2001). 

[5] A. Lotka, Proc. Natl. Acad. Sci. U.S.A 6, 410 (1920). 

[6] V. Volterra, Lecons sur la theorie mathematique de la 
lutte pour la vie (Gauthier-Villars, 1931). 

[7] T. Reichenbach, M. Mobilia, and E. Frey, Phys. Rev. E 
74, 051907 (2006). 

[8] M. Dykman, E. Mori, J. Ross, and P. Hunt, J. Chem. 
Phys. 100, 5735 (1994). 

[9] V. Elgart and A. Kamenev, Phys. Rev. E 70, 041106 
(2004). 

[10] A. Kamenev and B. Meerson, Phys. Rev. E 77, 061107 
(2008). 



[11] M. Dykman, I. Schwartz, and A. Landsman, Phys. Rev. 

Lett. 101, 078101 (2009). 
[12] B. Gaveau, M. Moreau, and J. Toth, Lett. Math. Phys. 

37, 285 (1996). 

[13] C. Doering, K. Sargsyan, and L. Sander, Multi-scale 

Model, and Simul. 3, 283 (2005). 
[14] M. Assaf and B. Meerson, Phys. Rev. Lett. 97, 200602 

(2007). 

[15] M. Assaf and B. Meerson, Phys. Rev. E 75, 031122 
(2007). 

[16] N. van Kampen, Stochastic Processes in Physics and 
Chemistry (North-Holland, Amsterdam, 2001). 

[17] D. Apalkov and P. Visscher, Phys. Rev. B 72, 180405 
(2005). 

[18] R. Durrett, SIAM Review 41, 677 (1999). 

[19] R. Abta, M. Schiffer, and N. Shnerb, Phys. Rev. Lett. 

98, 098104 (2007). 
[20] M. Mobilia, I. Georgiev, and U. Tauber, Jour. Stat. Phys. 

128, 447 (2007). 



